using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.IO;
using System.Text;
using System.Text.RegularExpressions;

using Atomic.Interfaces;
using Atomic.Interfaces.Atat;
using Atomic.Libraries;
using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Mathematics.RandomNumbers;
using Atomic.Structures.Symmetries;
using Atomic.Vasp;

namespace Atomic.Structures
{
	/// <summary>
	/// Structure class representing a unit cell of atoms. Immutable by design; manipulations are performed on copies of the structures.
	/// </summary>
	[Serializable]
	[DebuggerDisplay("{Formula}")]
	public class Structure
	{
		[NonSerialized] int hashCode;

		/// <summary>
		/// Creates a new structure using the lattice as the basis for the coordinate system.
		/// </summary>
		public Structure(Lattice lattice, params Site[] sites)
			: this(lattice, (IEnumerable<Site>)sites)
		{
		}

		/// <summary>
		/// Creates a new structure using the lattice as the basis for the coordinate system.
		/// </summary>
		public Structure(Lattice lattice, IEnumerable<Site> sites)
			: this(lattice, lattice, sites)
		{
		}

		/// <summary>
		/// Creates a new structure with explicitly specified basis for the coordinate system. If no basis is specified, physical coordinates are used.
		/// </summary>
		public Structure(Basis basis, Lattice lattice, params Site[] sites)
			: this(basis, lattice, (IEnumerable<Site>)sites)
		{
		}
		
		/// <summary>
		/// Creates a new structure with explicitly specified basis for the coordinate system. If no basis is specified, physical coordinates are used.
		/// </summary>
		public Structure(Basis basis, Lattice lattice, IEnumerable<Site> sites)
		{	
			if (basis == null || lattice == null || sites == null)
			{
				throw new ArgumentNullException();
			}

			Basis = basis;
			Lattice = lattice;
			ReciprocalLattice = lattice.ReciprocalLattice;
			Volume = lattice.Volume;
			Sites = new ReadOnlyList<Site>(sites);
			Formula = BuildFormula();
			SpaceGroup = new SpaceGroup(this);
		}

		/// <summary>
		/// Creates a new structure given atoms and relative positions.
		/// </summary>
		public Structure(Lattice lattice, IEnumerable<Atom> atoms, IEnumerable<SpaceVector> positions)
			: this(lattice, CreateCellAtoms(lattice, new List<Atom>(atoms).ToArray(), new List<SpaceVector>(positions).ToArray()))
		{
		}

		/// <summary>
		/// Creates a new structure given a formula (e.g. "H2O" for water) and relative positions of the atoms (in the same order as the appearance in the formula).
		/// </summary>
		public Structure(Lattice lattice, string formula, params SpaceVector[] positions)
			: this(lattice, ParseFormula(formula), positions)
		{
		}

		private static Atom[] ParseFormula(string formula)
		{
			List<Atom> atoms = new List<Atom>();

			while (formula.Length > 0)
			{
				Match m = Regex.Match(formula, "^([A-Z][a-z]*)([1-9][0-9]*)?");
				if (!m.Success)
				{
					throw new ArgumentException("Unrecognized or invalid chemical formula.");
				}

				string symbol = m.Groups[1].Value;
				int n = m.Groups[2].Success ? int.Parse(m.Groups[2].Value) : 1;

				for (int i = 0; i < n; i++)
				{
					atoms.Add(Atom.Parse(symbol));
				}

				formula = formula.Substring(m.Length);
			}

			return atoms.ToArray();
		}

		private static Site[] CreateCellAtoms(Basis basis, Atom[] atoms, SpaceVector[] positions)
		{
			if (positions == null)
			{
				throw new ArgumentNullException();
			}

			int n = atoms.Length;

			if (positions.Length != n)
			{
				throw new ArgumentException("Number of atoms and positions disagree.");
			}

			Site[] sites = new Site[n];
			for (int i = 0; i < n; i++)
			{
				sites[i] = new Site(atoms[i], new Position(basis, positions[i]));
			}

			return sites;
		}

		/// <summary>
		/// Reads structure in ATAT format.
		/// </summary>
		public static Structure Read(FileInfo file)
		{
			return AtatStructure.Read(file);
		}

		/// <summary>
		/// Reads structure in ATAT format.
		/// </summary>
		public static Structure Read(TextReader reader)
		{
			return AtatStructure.Read(reader);
		}

		/// <summary>
		/// Reads structure in ATAT format.
		/// </summary>
		public static Structure Read(string s)
		{
			return AtatStructure.Read(s);
		}

		/// <summary>
		/// Writes the structure in ATAT format.
		/// </summary>
		public void Write(FileInfo file)
		{
			AtatStructure.Write(this, file);
		}

		/// <summary>
		/// Writes the structure in ATAT format.
		/// </summary>
		public void Write(TextWriter writer)
		{
			AtatStructure.Write(this, writer);
		}

		/// <summary>
		/// Writes the structure in ATAT format.
		/// </summary>
		public string Write()
		{
			return AtatStructure.Write(this);
		}

		/// <summary>
		/// Shows the structures using VESTA.
		/// </summary>
		public static void Show(IEnumerable<Structure> structures)
		{
			Vesta.ShowStructure(structures);
		}

		/// <summary>
		/// Shows the structures using VESTA.
		/// </summary>
		public static void Show(params Structure[] structures)
		{
			Show((IEnumerable<Structure>)structures);
		}

		/// <summary>
		/// Shows the structure using VESTA.
		/// </summary>
		public void Show()
		{
			Structure.Show(this);
		}

		/// <summary>
		/// Returns a copy of the structure with an atom added.
		/// </summary>
		public Structure Add(Site site)
		{
			List<Site> sites = new List<Site>(Sites);
			sites.Add(site);

			return new Structure(Lattice, sites);
		}

		public Structure Add(Atom atom, Position position)
		{
			// This overload allows using any basis independent of the lattice (e.g. physical coordinates).
			return Add(new Site(atom, position));
		}

		public Structure Add(Atom atom, SpaceVector coordinates)
		{
			return Add(atom, new Position(Lattice, coordinates));
		}

		public Structure Add(Atom atom, double x, double y, double z)
		{
			return Add(atom, new Position(Lattice, x, y, z));
		}
		
		/// <summary>
		/// Returns a copy of the structure with sites satisfying the predicate removed.
		/// </summary>
		public Structure Remove(Predicate<Site> match)
		{
			return new Structure(Lattice, new List<Site>(Sites).FindAll(site => !match(site)));
		}

		/// <summary>
		/// Returns a copy of the structure with specified atom removed.
		/// </summary>
		public Structure Remove(Atom atom)
		{
			return Remove(site => site.Atom == atom).Transform(a => a is PlaceholderAtom ? ((PlaceholderAtom)a).Remove(atom) : a);
		}

		/// <summary>
		/// Returns a copy of the structure with sites containing <see cref="VacancyAtom" /> removed.
		/// </summary>
		public Structure RemoveVacancies()
		{
			return Remove(site => site.Atom is VacancyAtom);
		}

		public Structure Transform(Func<Site, Site> transform)
		{
			return new Structure(Lattice, new List<Site>(Sites).ConvertAll<Site>(site => transform(site)));
		}

		public Structure Transform(Func<Atom, Atom> transform)
		{
			return Transform(site => new Site(transform(site.Atom), site.Position));
		}
		
		public Structure Transform(Func<Position, Position> transform)
		{
			return Transform(site => new Site(site.Atom, transform(site.Position)));
		}

		public Structure TransformShifted(SpaceVector coordinates)
		{
			List<Site> sites = new List<Site>();
			foreach (Site site in Sites)
			{
				// Perform shift in basis coordinates, but transform to lattice coordinates after this.
				SpaceVector p = Lattice.DecomposeRounded(new Position(Basis, Basis.Decompose(site.Position).Coordinates + coordinates)).Coordinates;

				// Subtracting an integral multiple of lattice vectors doesn't change the position.
				double x = p[0] - Math.Floor(p[0]);
				double y = p[1] - Math.Floor(p[1]);
				double z = p[2] - Math.Floor(p[2]);

				/*// Allow numbers very close to 1 being represented as 0 instead.
				x = Math.Abs(x - 1.0) < 1.0e-12 ? 0.0 : x;
				y = Math.Abs(y - 1.0) < 1.0e-12 ? 0.0 : y;
				z = Math.Abs(z - 1.0) < 1.0e-12 ? 0.0 : z;*/
				p = new SpaceVector(x, y, z);

				// Transform back to basis coordinates and add site.
				sites.Add(new Site(site.Atom, new Position(Basis, Basis.Decompose(new Position(Lattice, p)).Coordinates)));
			}

			return new Structure(Basis, Lattice, sites);
		}

		public Structure TransformShifted(double x, double y, double z)
		{
			return TransformShifted(new SpaceVector(x, y, z));
		}

		/// <summary>
		/// Creates a superstructure by repeating the structure an integer times along each lattice vector.
		/// </summary>
		public Structure Repeat(int m1, int m2, int m3)
		{
			if (m1 < 1 || m2 < 1 || m3 < 1)
			{
				throw new ArgumentException("Can only repeat a cell by a positive integer multiple.");
			}

			//return Repeat(m1, 0, 0, 0, m2, 0, 0, 0, m3);

			// Keep the parent basis unchanged (if defined) while enlarging the lattice vectors.
			SpaceVector a = Basis.Decompose(Lattice.A).Coordinates;
			SpaceVector b = Basis.Decompose(Lattice.B).Coordinates;
			SpaceVector c = Basis.Decompose(Lattice.C).Coordinates;
			Lattice lattice = Lattice.Parse(Basis, m1 * a, m2 * b, m3 * c);

			List<Site> sites = new List<Site>();
			for (int i = 0; i < m1; i++)
			{
				for (int j = 0; j < m2; j++)
				{
					for (int k = 0; k < m3; k++)
					{
						foreach (Site site in Sites)
						{
							// Make sure the coordinates are relative to the lattice of the original structure.
							SpaceVector coordinates = Basis.Decompose(site.Position).Coordinates;
							sites.Add(new Site(site.Atom, new Position(Basis, coordinates + i * a + j * b + k * c)));
						}
					}
				}
			}

			return new Structure(Basis, lattice, sites);
		}

		/// <summary>
		/// Creates a superstructure by repeating the structure an integer times along each lattice vector.
		/// </summary>
		public Structure Repeat(int m11, int m12, int m13, int m21, int m22, int m23, int m31, int m32, int m33)
		{
			if (m11 >= 1 && m22 >= 1 && m33 >= 1 && m12 == 0 && m13 == 0 && m21 == 0 && m23 == 0 && m31 == 0 && m32 == 0)
			{
				// Handle this case with the closed form solution above.
				return Repeat(m11, m22, m33);
			}

			// Keep the parent basis unchanged (if defined) while enlarging the lattice vectors.
			SpaceVector a = Basis.Decompose(Lattice.A).Coordinates;
			SpaceVector b = Basis.Decompose(Lattice.B).Coordinates;
			SpaceVector c = Basis.Decompose(Lattice.C).Coordinates;

			// The lattice vectors of the supercell is defined as column vectors in the relation
			// (as,bs,cs)=(a,b,c)*M, like described here:
			// http://phonopy.sourceforge.net/setting-tags.html#dim
			//SpaceMatrix m = new SpaceMatrix(m11, m12, m13, m21, m22, m23, m31, m32, m33);

			// Linear independence of lattice vectors is enforced by the Lattice class constructor.
			Lattice lattice = Lattice.Parse(Basis, m11 * a + m21 * b + m31 * c, m12 * a + m22 * b + m32 * c, m13 * a + m23 * b + m33 * c);

			// Implement a better method than this naive brute force method!!
			List<Site> sites = new List<Site>();
			for (int i = -10; i < 10; i++)
			{
				for (int j = -10; j < 10; j++)
				{
					for (int k = -10; k < 10; k++)
					{
						foreach (Site site in Sites)
						{
							// Make sure the coordinates are relative to the lattice of the original structure.
							SpaceVector coordinates = Basis.Decompose(site.Position).Coordinates;

							// Include only if within the standard unit cell.
							SpaceVector p = lattice.DecomposeRounded(new Position(Basis, coordinates + i * a + j * b + k * c)).Coordinates;

							double x = p[0];
							double y = p[1];
							double z = p[2];

							/*// Rounding around 0.
							x = Math.Abs(x) < 1.0e-12 ? 0.0 : x;
							y = Math.Abs(y) < 1.0e-12 ? 0.0 : y;
							z = Math.Abs(z) < 1.0e-12 ? 0.0 : z;

							// Rounding around 1.
							x = Math.Abs(x - 1.0) < 1.0e-12 ? 1.0 : x;
							y = Math.Abs(y - 1.0) < 1.0e-12 ? 1.0 : y;
							z = Math.Abs(z - 1.0) < 1.0e-12 ? 1.0 : z;*/

							//if (x + 1.0e-10 >= 0.0 && x + 1.0e-10 < 1.0 && y + 1.0e-10 >= 0.0 && y + 1.0e-10 < 1.0 && z + 1.0e-10  >= 0.0 && z + 1.0e-10 < 1.0)
							if (x >= 0.0 && x < 1.0 && y >= 0.0 && y < 1.0 && z >= 0.0 && z < 1.0)
							{
								sites.Add(new Site(site.Atom, new Position(lattice, new SpaceVector(x, y, z))));
							}
						}
					}
				}
			}

			// Check expected scaling.
			int n = (int)Math.Round(lattice.Volume / Lattice.Volume);
			if (sites.Count != n * Sites.Count)
			{
				throw new Exception("Couldn't build repeated structure.");
			}

			return new Structure(Basis, lattice, sites);
		}

		/// <summary>
		/// Creates a superstructure by repeating the structure an integer times along each lattice vector.
		/// </summary>
		public Structure Repeat(LatticeRepeatMatrix m)
		{
			return Repeat(m[0, 0], m[0, 1], m[0, 2], m[1, 0], m[1, 1], m[1, 2], m[2, 0], m[2, 1], m[2, 2]);
		}

		/// <summary>
		/// Creates a superstructure by specifying the minimum distance of between mirrored images of the same original atom.
		/// </summary>
		public Structure Repeat(double distance)
		{
			return AtatSuperStructureGenerator.GenerateMapped(this, distance);
		}

		/// <summary>
		/// Tries to rebuild the original structure from a superstructure.
		/// </summary>
		public Structure InverseRepeat(int m1, int m2, int m3)
		{
			if (m1 < 1 || m2 < 1 || m3 < 1)
			{
				throw new ArgumentException("Can only repeat a cell by a positive integer multiple.");
			}

			// Keep the parent basis unchanged (if defined) while enlarging the lattice vectors.
			SpaceVector a = Basis.Decompose(Lattice.A).Coordinates;
			SpaceVector b = Basis.Decompose(Lattice.B).Coordinates;
			SpaceVector c = Basis.Decompose(Lattice.C).Coordinates;
			Lattice lattice = Lattice.Parse(Basis, a / m1, b / m2, c / m3);

			// Implement a better method than this naive brute force method!!
			List<Site> sites = new List<Site>();

			foreach (Site site in Sites)
			{
				// Make sure the coordinates are relative to the lattice of the original structure.
				SpaceVector coordinates = Basis.Decompose(site.Position).Coordinates;

				// Include only if within the standard unit cell.
				SpaceVector p = lattice.DecomposeRounded(new Position(Basis, coordinates)).Coordinates;

				double x = p[0];
				double y = p[1];
				double z = p[2];

				if (x >= 0.0 && x < 1.0 && y >= 0.0 && y < 1.0 && z >= 0.0 && z < 1.0)
				{
					sites.Add(new Site(site.Atom, new Position(lattice, new SpaceVector(x, y, z))));
				}
			}

			// Check expected scaling.
			//int n = (int)Math.Round(Lattice.Volume / lattice.Volume);
			int n = m1 * m2 * m3;
			if (n * sites.Count != Sites.Count)
			{
				throw new Exception("Couldn't determine inverse of superstructure.");
			}

			return new Structure(Basis, lattice, sites);
		}

		/// <summary>
		/// Returns a copy of the structure scaled uniformly in each direction.
		/// </summary>
		public Structure Scale(double scaling)
		{
			// Scale basis, but keep relative lattice vectors. Allow for the case where the original basis is undefined.
			SpaceVector a = Basis.Decompose(Lattice.A).Coordinates;
			SpaceVector b = Basis.Decompose(Lattice.B).Coordinates;
			SpaceVector c = Basis.Decompose(Lattice.C).Coordinates;
			Basis scaledBasis = Basis.Scale(scaling);

			// Define lattice vectors relative to this scaled basis.
			Lattice lattice = Lattice.Parse(scaledBasis, a, b, c);

			List<Site> sites = new List<Site>();
			foreach (Site site in Sites)
			{
				// Make sure the coordinates are relative to the lattice of the original structure. Then keep
				// the relative coordinates, but change the lattice.
				SpaceVector coordinates = Basis.Decompose(site.Position).Coordinates;
				sites.Add(new Site(site.Atom, new Position(scaledBasis, coordinates)));
			}
	
			return new Structure(scaledBasis, lattice, sites);
		}

		/// <summary>
		/// Returns a copy of the structure with the lattice vectors scaled.
		/// </summary>
		public Structure Scale(double s1, double s2, double s3)
		{
			if (s1 <= 0.0 || s2 <= 0.0 || s3 <= 0.0)
			{
				throw new ArgumentOutOfRangeException();
			}

			SpaceVector a = s1 * Lattice.A.PhysicalCoordinates;
			SpaceVector b = s2 * Lattice.B.PhysicalCoordinates;
			SpaceVector c = s3 * Lattice.C.PhysicalCoordinates;
			Lattice latticeScaled = new Lattice(Basis.StandardBasis, a, b, c);

			List<Site> sites = new List<Site>();
			foreach (Site site in Sites)
			{
				// Make sure the coordinates are relative to the lattice of the original structure. Then keep
				// the relative coordinates, but change the lattice.
				SpaceVector coordinates = Lattice.Decompose(site.Position).Coordinates;
				sites.Add(new Site(site.Atom, new Position(latticeScaled, coordinates)));
			}

			return new Structure(latticeScaled, sites);
		}

		/// <summary>
		/// Returns a copy of the structure with the volume scaled.
		/// </summary>
		public Structure ScaleVolume(double scaling)
		{
			return Scale(Math.Pow(scaling, 1.0 / 3.0));
		}

		public Structure Rotate(SpaceMatrix rotation)
		{
			// Rotate basis vectors and keep lattice vectors defined relative to the basis.
			SpaceVector a = Basis.Decompose(Lattice.A).Coordinates;
			SpaceVector b = Basis.Decompose(Lattice.B).Coordinates;
			SpaceVector c = Basis.Decompose(Lattice.C).Coordinates;
			Basis rotatedBasis = Basis.Rotate(rotation);

			// Define lattice vectors relative to this scaled basis.
			Lattice lattice = Lattice.Parse(rotatedBasis, a, b, c);

			List<Site> sites = new List<Site>();
			foreach (Site site in Sites)
			{
				// Make sure the coordinates are relative to the lattice of the original structure. Then keep
				// the relative coordinates, but change the lattice.
				SpaceVector coordinates = Basis.Decompose(site.Position).Coordinates;
				sites.Add(new Site(site.Atom, new Position(rotatedBasis, coordinates)));
			}
	
			return new Structure(rotatedBasis, lattice, sites);
		}

		public Structure SpinUnpolarize()
		{
			return Transform(atom => atom.SpinUnpolarize());
		}

		public Structure Normalize()
		{
			List<Site> sites = new List<Site>();
			foreach (Site site in Sites)
			{
				sites.Add(new Site(site.Atom, Lattice.DecomposeRounded(site.Position)));
			}

			return new Structure(Basis, Lattice, sites);
		}

		/// <summary>
		/// Copies the structure with atoms from the first placeholders being replaced by the specified atoms.
		/// </summary>
		public Structure FillPlaceholders(IEnumerable<Atom> atoms)
		{
			return FillPlaceholders(new List<Atom>(atoms).ToArray());
		}

		/// <summary>
		/// Copies the structure with atoms from the first placeholders being replaced by the specified atoms.
		/// </summary>
		public Structure FillPlaceholders(params Atom[] atoms)
		{
			int i = 0;

			List<Site> sites = new List<Site>();
			foreach (Site site in Sites)
			{
				if (i < atoms.Length && site.Atom is PlaceholderAtom)
				{
					// Allow the placeholder to be kept. Otherwise the placeholder has to contain the specified atom.
					if (site.Atom != atoms[i] && !((PlaceholderAtom)site.Atom).Atoms.Contains(atoms[i]))
					{
						throw new Exception("Placeholder doesn't contain specified atom.");
					}

					sites.Add(new Site(atoms[i], site.Position));
					i++;
				}
				else
				{
					sites.Add(site);
				}
			}

			if (i != atoms.Length)
			{
				throw new Exception("More atoms specified than available placeholders.");
			}

			return new Structure(Basis, Lattice, sites);
		}
		
		/// <summary>
		/// Copies the structure with atoms from the first placeholders being replaced by the specified atom.
		/// </summary>
		public Structure FillPlaceholders(Atom atom, int count)
		{
			if (count < 0)
			{
				throw new ArgumentException();
			}

			Atom[] atoms = new Atom[count];
			for (int i = 0; i < count; i++)
			{
				atoms[i] = atom;
			}

			return FillPlaceholders(atoms);
		}

		/// <summary>
		/// Copies the structure with atoms from placeholders being selected randomly.
		/// </summary>
		public Structure FillPlaceholders(IRandomGenerator randomGenerator)
		{
			return Transform(atom => atom is PlaceholderAtom ? ((PlaceholderAtom)atom).Atoms[randomGenerator.Next(((PlaceholderAtom)atom).Atoms.Count)] : atom);
		}

		/// <summary>
		/// Fills placeholders with the corresponding atoms in all possible ways, not considering symmetry. The may be faster than trying to figure out symmetrically district combinations.
		/// </summary>
		public IEnumerable<Structure> FillPlaceholdersCombinatorial()
		{
			return StructureGeneratorNaive.Generate(this);
		}

		/// <summary>
		/// Fills placeholders with the corresponding atoms in all possible symmetrically distinct ways. A list of resulting structures is returned.
		/// This uses the ATAT genstr tool.
		/// </summary>
		public IEnumerable<Structure> FillPlaceholdersCombinatorialReduced()
		{
			return AtatStructureGenerator.GenerateMapped(this);
		}

		/// <summary>
		/// Fills placeholders a number of placeholder with the specified atom in all possible ways, not considering symmetry. The remaining placeholders are left unchanged.
		/// </summary>
		public IEnumerable<Structure> FillPlaceholdersCombinatorialPartial(Atom atom, int count)
		{
			return StructureGeneratorNaive.GeneratePartial(this, atom, count);
		}

		/// <summary>
		/// Builds the chemical formula for the atoms in the structure. Repeated atoms are automatically grouped and counted. A preferred order may optionally be specified.
		/// </summary>
		public string BuildFormula(params Atom[] order)
		{
			List<Atom> atoms = new List<Atom>(order);
			Dictionary<Atom, int> counts = new Dictionary<Atom, int>();

			foreach (Site site in Sites)
			{
				if (site == null)
				{
					continue;
				}

				Atom atom = site.Atom;
				if (!atoms.Contains(atom))
				{
					atoms.Add(atom);
				}
				if (!counts.ContainsKey(atom))
				{
					counts.Add(atom, 0);
				}
				counts[atom] = counts[atom] + 1;
			}

			StringBuilder formula = new StringBuilder();
			foreach (Atom atom in atoms)
			{
				if (counts.ContainsKey(atom))
				{
					string symbol = atom.Symbol;
					if (atom is PlaceholderAtom || atom is SpinPolarizedAtom)
					{
						symbol = "(" + symbol + ")";
					}
					if (counts[atom] > 1)
					{
						symbol += counts[atom];
					}
					formula.Append(symbol);

					// Skip this atom next time.
					counts.Remove(atom);
				}
			}

			return formula.ToString();
		}

		public virtual bool Equals(Structure other)
		{
			if (object.ReferenceEquals(other, null))
			{
				return false;
			}

			if (other is VaspStructure)
			{
				// Special rules here. Implemented in the other class.
				return other.Equals(this);
			}

			return EqualityHelper.Equals(this, other, new object[] { Lattice, new List<Site>(Sites).ToArray() }, new object[] { other.Lattice, new List<Site>(other.Sites).ToArray() }, hashCode, other.hashCode);
		}

		public override bool Equals(object other)
		{
			return Equals(other as Structure);
		}
		
		public override int GetHashCode()
		{
			if (hashCode == 0)
			{
				hashCode = EqualityHelper.HashCode(new object[] { Lattice, new List<Site>(Sites).ToArray() });
			}

			return hashCode;
		}

		/// <summary>
		/// Compares two structures allowing small deviation in the atomic positions.
		/// This does NOT define an equivalence relation as it lacks the transitive property.
		/// The tolerance is assumed to be in real units (i.e. in Å).
		/// </summary>
		public bool Equals(Structure other, double tolerance)
		{
			return Equals(this, other, tolerance);
		}

		/// <summary>
		/// Compares two structures allowing small deviation in the atomic positions.
		/// This does NOT define an equivalence relation as it lacks the transitive property.
		/// The tolerance is assumed to be in real units (i.e. in Å).
		/// </summary>
		public static bool Equals(Structure structure1, Structure structure2, double tolerance)
		{
			return StructureEquality.Equals(structure1, structure2, tolerance);
		}

		/// <summary>
		/// Compares two structures allowing small deviation in the atomic positions and acting with symmetry operations.
		/// This does NOT define an equivalence relation as it lacks the transitive property.
		/// The tolerance is assumed to be in real units (i.e. in Å).
		/// </summary>
		public bool SymmetryEquals(Structure other, IEnumerable<SpaceGroupSymmetryOperation> symmetries, double tolerance)
		{
			return SymmetryEquals(this, other, symmetries, tolerance);
		}

		/// <summary>
		/// Compares two structures allowing small deviation in the atomic positions and acting with symmetry operations.
		/// This does NOT define an equivalence relation as it lacks the transitive property.
		/// The tolerance is assumed to be in real units (i.e. in Å).
		/// </summary>
		public static bool SymmetryEquals(Structure structure1, Structure structure2, IEnumerable<SpaceGroupSymmetryOperation> symmetries, double tolerance)
		{
			return StructureEquality.SymmetryEquals(structure1, structure2, symmetries, tolerance);
		}

		/// <summary>
		/// Computes the maximum distance between sites in the two structures. Lattice vectors are allowed to be different.
		/// </summary>
		public double Distance(Structure other)
		{
			return Distance(this, other);
		}

		/// <summary>
		/// Computes the maximum distance between sites in the two structures. Lattice vectors are allowed to be different.
		/// </summary>
		public static double Distance(Structure structure1, Structure structure2)
		{
			return StructureEquality.Distance(structure1, structure2);
		}

		/// <summary>
		/// Creates an empty simple cubic structure.
		/// </summary>
		public static Structure CreateSimpleCubicStructure(double a)
		{
			return CreateSimpleCubicStructure(a, new Atom[0], new SpaceVector[0]);
		}

		/// <summary>
		/// Creates a simple cubic structure.
		/// </summary>
		public static Structure CreateSimpleCubicStructure(double a, IEnumerable<Atom> atoms, IEnumerable<SpaceVector> positions)
		{
			//return new Structure(new SimpleCubicLattice(a), atoms, positions);
			return new Structure(Lattice.CreateSimpleCubicLattice(a), atoms, positions);
		}

		/// <summary>
		/// Creates a simple cubic structure.
		/// </summary>
		public static Structure CreateSimpleCubicStructure(double a, string formula, params SpaceVector[] positions)
		{
			return CreateSimpleCubicStructure(a, ParseFormula(formula), positions);
		}

		/// <summary>
		/// Creates an empty fcc structure.
		/// </summary>
		public static Structure CreateFaceCenteredCubicStructure(double a)
		{
			return CreateFaceCenteredCubicStructure(a, new Atom[0], new SpaceVector[0]);
		}

		/// <summary>
		/// Creates a fcc structure.
		/// </summary>
		public static Structure CreateFaceCenteredCubicStructure(double a, IEnumerable<Atom> atoms, IEnumerable<SpaceVector> positions)
		{
			return new Structure(new FaceCenteredCubicLattice(a), atoms, positions);
		}

		/// <summary>
		/// Creates a fcc structure.
		/// </summary>
		public static Structure CreateFaceCenteredCubicStructure(double a, string formula, params SpaceVector[] positions)
		{
			return CreateFaceCenteredCubicStructure(a, ParseFormula(formula), positions);
		}

		/// <summary>
		/// Creates an empty bcc structure.
		/// </summary>
		public static Structure CreateBodyCenteredCubicStructure(double a)
		{
			return CreateBodyCenteredCubicStructure(a, new Atom[0], new SpaceVector[0]);
		}

		/// <summary>
		/// Creates a bcc structure.
		/// </summary>
		public static Structure CreateBodyCenteredCubicStructure(double a, IEnumerable<Atom> atoms, IEnumerable<SpaceVector> positions)
		{
			return new Structure(new BodyCenteredCubicLattice(a), atoms, positions);
		}

		/// <summary>
		/// Creates a bcc structure.
		/// </summary>
		public static Structure CreateBodyCenteredCubicStructure(double a, string formula, params SpaceVector[] positions)
		{
			return CreateBodyCenteredCubicStructure(a, ParseFormula(formula), positions);
		}

		/// <summary>
		/// Creates an empty hcp structure.
		/// </summary>
		public static Structure CreateHexagonalStructure(double a)
		{
			return CreateHexagonalStructure(a, new Atom[0], new SpaceVector[0]);
		}

		/// <summary>
		/// Creates a hcp structure.
		/// </summary>
		public static Structure CreateHexagonalStructure(double a, IEnumerable<Atom> atoms, IEnumerable<SpaceVector> positions)
		{
			// Ideal ratio 3c^2=8a^2.
			return CreateHexagonalStructure(a, a * 1.6329931618554521, atoms, positions);
		}

		/// <summary>
		/// Creates a hcp structure.
		/// </summary>
		public static Structure CreateHexagonalStructure(double a, string formula, params SpaceVector[] positions)
		{
			return CreateHexagonalStructure(a, ParseFormula(formula), positions);
		}

		/// <summary>
		/// Creates an empty hcp structure with ideal lattice plane ratio.
		/// </summary>
		public static Structure CreateHexagonalStructure(double a, double c)
		{
			return CreateHexagonalStructure(a, c, new Atom[0], new SpaceVector[0]);
		}

		/// <summary>
		/// Creates a hcp structure with ideal lattice plane ratio.
		/// </summary>
		public static Structure CreateHexagonalStructure(double a, double c, IEnumerable<Atom> atoms, IEnumerable<SpaceVector> positions)
		{
			//return new Structure(new HexagonalLattice(a, c), atoms, positions);
			return new Structure(Lattice.CreateHexagonalLattice(a, c), atoms, positions);
		}

		/// <summary>
		/// Creates a hcp structure with ideal lattice plane ratio.
		/// </summary>
		public static Structure CreateHexagonalStructure(double a, double c, string formula, params SpaceVector[] positions)
		{
			return CreateHexagonalStructure(a, c, ParseFormula(formula), positions);
		}

		public static bool operator ==(Structure structure1, Structure structure2)
		{
			if (object.ReferenceEquals(structure1, null))
			{
				if (object.ReferenceEquals(structure2, null))
				{
					return true;
				}

				return false;
			}

			return structure1.Equals(structure2);
		}

		public static bool operator !=(Structure structure1, Structure structure2)
		{
			return !(structure1 == structure2);
		}

		/// <summary>
		/// Default basis for coordinates when adding new sites.
		/// </summary>
		public Basis Basis
		{
			get;
			private set;
		}

		/// <summary>
		/// Lattice vectors defining the translational periodicity of the structure.
		/// </summary>
		public Lattice Lattice
		{
			get;
			private set;
		}

		/// <summary>
		/// Reciprocal lattice vectors.
		/// </summary>
		public Lattice ReciprocalLattice
		{
			get;
			private set;
		}

		/// <summary>
		/// Unit cell volume.
		/// </summary>
		public double Volume
		{
			get;
			private set;
		}

		/// <summary>
		/// The atomic sites (i.e. positions and species of atoms) of the structure.
		/// </summary>
		public ReadOnlyList<Site> Sites
		{
			get;
			private set;
		}

		/// <summary>
		/// The chemical formula for the atoms in the structure.
		/// </summary>
		public string Formula
		{
			get;
			private set;
		}

		public SpaceGroup SpaceGroup
		{
			get;
			private set;
		}
	}
}
